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Abstract 



Estimation of the operational risk capital under the Loss Distribution 
Approach requires evaluation of aggregate (compound) loss distribu- 
tions which is one of the classic problems in risk theory. Closed-form 
solutions arc not available for the distributions typically used in oper- 
ational risk. However with modern computer processing power, these 
^ distributions can be calculated virtually exactly using numerical meth- 

00 ods. This paper reviews numerical algorithms that can be successfully 

used to calculate the aggregate loss distributions. In particular Monte 
I Carlo, Panjer recursion and Fourier transformation methods are pre- 

^ sented and compared. Also, several closed-form approximations based 

• ^ on moment matching and asymptotic result for heavy-tailed distribu- 

^ tions are reviewed. 



Keywords: aggregate loss distribution, compound distribution, Monte 
Carlo, Panjer recursion, Fast Fourier Transform, loss distribution ap- 
proach, operational risk. 



1 Introduction and Model 

Estimation of the operational risk capital under the Loss Distribution Ap- 
proach (LDA) requires calculation of the distribution for the aggregate (com- 
pound) loss 

Z^X, + --- + Xn, (1) 
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where the frequency is a discrete random variable and Xi, . . . jX^r are 
positive random severities. For a recent review of LDA, see Chernobai et 
al (2007) and Shevchenko (2010). This is one of the classical problems in 
risk theory. Closed-form solutions are not available for the distributions typ- 
ically used in operational risk. However with modern computer processing 
power, these distributions can be calculated virtually exactly using numerical 
algorithms. The easiest to implement is the Monte Carlo method. However, 
because it is typically slow, Panjer recursion and Fourier inversion techniques 
are widely used. Both have a long history, but their applications to comput- 
ing very high quantiles of the compound distribution functions with high 
frequencies and heavy tails are only recent developments and various pitfalls 
exist. 

This paper presents review and tutorial on the methods used to calculate 
the distribution of the aggregate loss ([T]) over a chosen time period. The 
following model assumptions and notation are used: 

• Only one risk cell and one time period are considered. Typically, the 
calculation of the aggregate loss over a one-year time period is required 
in operational risk. 

• N is the number of events over the time period (frequency) modelled 
as a discrete random variable with probability mass function pk = 
Pi[N = k], k = 0,1,... . There is a finite probability of no loss 
occurring over the considered time period if X = is allowed, i.e. 
Pt[Z = 0] = Fr[N = 0]. 

• Xj, i > 1 are positive severities of the events (loss amounts) modelled 
as independent and identically distributed random variables from a 
continuous distribution function F{x) with a; > and F{0) = 0. The 
corresponding density function is denoted as f{x). 

• N and Xi are independent for all i, i.e. the frequencies and severities 
are independent. 

• The distribution and density functions of the aggregate loss Z are de- 
noted as H{z) and h{z) respectively. 

• All model parameters (parameters of the frequency and severity dis- 
tributions) are assumed to be known. In real application, the model 
parameters are unknown and estimated using past data. The impact 
of uncertainty in parameter estimates on the annual loss distribution 
can be significant for low-frequency/high-severity operational risks due 
to limited historical data (see Shevchenko (2008)); this topic is beyond 
the purpose of this paper. 

In general, there are two types of analytic solutions for calculating the 
compound distribution H{z). These are based on convolutions and method 
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of characteristic functions described in Section IH The moments of the com- 
pound loss can be derived in closed-form via the moments of frequency and 
severity; these are presented in Section [2] as well. Section [3] gives the analytic 
expressions for the Value-at-Risk and expected shortfall risk measures. Typ- 
ically, the analytic solutions do not have closed-form and numerical methods 
such as Monte Carlo (MC), Panjer recursion, Fast Fourier Transform (FFT) 
or direct integration are required; these are described in Sections|4}[5|[6]and[7] 
respectively. Comparison of these methods is discussed in Section |8| Finally, 
Section |9] reviews several closed-form approximations. The distributions used 
throughout the paper are formally defined in Appendix. 



2 Analytic Solutions 

Analytic calculation of the compound distribution can be accomplished using 
methods of convolutions and characteristic functions. This section presents 
these methods and derives the moments of the compound distribution. 



2.1 Solution via Convolutions 

It is well-known that the density and distribution functions of the sum of two 
independent continuous random variables Yi ~ and I2 ~ -^2(")) with 

the densities /i(-) and /2(-) respectively, can be calculated via convolution 

as 

fn+vM = ifi * f2){y) = J f2{y - yi)fi{yi)dyi (2) 

and 

Fy,-,yM = (Fi * F2){y) = J F^iy - yi)Myi)dy, (3) 

respectively. Hereafter, notation fi * /2 denotes convolution of /i and /2 
functions as defined above; notation Y ~ F{y) means a random variable Y 
has a distribution function F{y). Thus the distribution of the aggregate loss 
([1]) can be calculated via convolutions as 



H{z) = Pt[Z < z] = '^Pt[Z < z\N = k]PT[N = k] 

00 

= (4) 

k=0 

Here, = Pr[Xi + ■ ■ ■ + Xk < z] is the k-th convolution of F(-) 

calculated recursively as 

F^^>{z)= r F^^-^>{z-x)f{x)dx 
Jo 

with 



_p(o)*| 



z 



1, z>0, 
0, z<0. 
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Note that the integration hmits are and z because the considered severities 
are nonnegative. Though the obtained formula is analytic, its direct calcula- 
tion is difficult because, in general, the convolution powers are not available 
in closed-form. Panjer recursion and FFT, discussed in Sections [5] and |6| are 
very efficient numerical methods to calculate these convolutions. 



2.2 Solution via Characteristic Functions 

The method of characteristic functions for computing probability distribu- 
tions is a powerful tool in mathematical finance; it is explained in many 
textbooks on probability theory. In particular, it is used for calculating ag- 
gregate loss distributions in the insurance, operational risk and credit risk. 
Typically, compound distributions cannot be found in closed-form but can 
be conveniently expressed through the inverse transform of the characteristic 
functions. The characteristic function of the severity density j{x) is formally 
defined as 

oo 

ip{t) = [ f{x)e''^dx, (5) 



where i = is a unit imaginary number. Also, the probability generating 
function of a frequency distribution with probability mass function = 
Ft[N = k] is 

oo 

Hs) = J2^%- (6) 

k=0 

Then, the characteristic function of the compound loss Z in model ([T]), de- 
noted by x{t)j can be expressed through the probability generating function 
of the frequency distribution and characteristic function of the severity dis- 
tribution as 

oo 
k=0 

For example: 

• If frequency is distributed from Poisson{X), then 

oo —A \ k 

X{t) = ^ = ^MMt) - A); (8) 

fe=0 

• If is from negative binomial distribution NegBin{m,p), then 

it) = f:(^(t))^Y'+7"M(i-P)v- 

I— n V / 



A;=0 

P 



l-{l-p)^{t) 



(9) 
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Given characteristic function, the density of the aggregate loss Z can be 
calculated via the inverse Fourier transform as 



oo 

h{z) = — x{t) exp{-itz)dt, z>0. (10) 
2n J 

— oo 

In the case of nonnegative severities, the density and distribution functions of 
the compound loss can be calculated using the following lemma (for a proof, 
see e.g. Luo and Shevchenko (2009, Appendix A)). 

Lemma 2.1 For a nonnegative random variable Z with a characteristic func- 
tion x(t), the density h{z) and distribution H{z) functions, z>{], are 



oo 



h{z) = - j Re[x(t)] cos{tz)dt, z>0; (11) 



oo 



H{z) = - [ Re[x{t)f-^^^^dt, z>0. (12) 



TT 



Changing variable x = t x z, the formula (12) can be rewritten as 



oo 



N 2 /"^ r / / Msin(x) , 
Hiz) = - j Re[xix/ z)]^^dx, 







which is often a useful representation to study hmiting properties. In partic- 
ular, in the limit z — )■ 0, it gives 



oo 



H{z ^ 0) = -Re[x(oo)] [ ^-^^dx = Re[x(oo)]. 

71 J X 



This leads to a correct limit H{0) = Pt[N = 0], because the severity char- 
acteristic function (p{oo) — 0. For example, H{0) = exp(— A) in the case of 
~ Poisson{X), and H{0) = for A^ ~ NegBin{m,p). 
FFT and direct integration methods to calculate the above Fourier trans- 
forms are discussed in details in Sections [6] and [7] respectively. 



2.3 Compound Distribution Moments 

In general, the compound distribution cannot be found in closed-form. How- 
ever, its moments can be expressed through the moments of the frequency 
and severity. It is convenient to calculate the moments via characteristic 
function. In particular, one can calculate the moments as 

, A; = l,2,.... (13) 

t=o 



dt^ 
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Similarly, the central moments can be found as 

/i, = E[(Z-E[Z])^] 

. d'x{t)exp{-ttE[Z]) 



, k = l,2,.... (14) 

i=0 



Here, for compound distribution, x(t) is given by ([T]). Then, one can derive 
the explicit expressions for all moments of compound distribution via the 
moments of frequency and severity noting that ip{0) = 1 and using relations 

^ = E[iV(iV-l)---(iV-fc + l)], (15) 
-ir^^ = E[Xf], (16) 



that follow from the definitions of the probability generating and character- 
istic functions ^ and ([s]) respectively, though the expression is lengthy for 
high moments. Sometimes, it is easier to work with the so-called cumulants 
(or semi-invariants) 



d^ Inx(t) 



(17) 

t=o 



dt^ 

which are closely related to the moments. The moments can be calculated 
via the cumulants and vice versa. In application, only the first four moments 
are most often used with the following relations: 

^2 = K2= Var[Z]; /is = k^; /i4 = K4 Snj. (18) 

Also, popular distribution characteristics are skewness = jJ^z/ {^12)^^'^ ^^"i 
kurtosis = — 3 + ii^/ {fi2Y- 

The above formulas relating characteristic function and moments can be 
found in many textbooks on risk theory such as McNeil et al (2005, Section 
10.2.2). The explicit expressions for the first four moments are given by the 
following proposition. 

Proposition 2.2 (Moments of compound distribution) The first four 
moments of the compound random variable Z = Xi + • ■ ■ -|- X^, where 
Xi,. . . ,Xn are independent and identically distributed, and independent of 
N , are given by 

E[Z] = E[X]E[Xi], 
Var[Z] = E[X]Var[Xi] + Var[X](E[Xi])2, 
E[{Z - E[Z]f] = E[X]E[(Xi - E[Xi]f] + 3Var[X]Var[Xi]E[Xi] 

+E[(X-E[X])3](E[Xl])^ 
E[(Z - E[Z])^] = E[X]E[(Xi - E[Xi])^] + 4Var[X]E[(Xi - E[Xi])3]E[Xi] 
+3(Var[X] + E[X](E[X] - l))(Var[Xi])2 
+6(E[(X - E[N]f] + E[X]Var[X])(E[Xi])2Var[Xi] 
+E[(X-E[X])^](E[Xi])^ 

Here, it is assumed that the required moments of severity and frequency exist. 
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Proof 1 This follows from the expression for characteristic function of the 
compound distribution ^ and formulas ( 15p6 ). The calculus is simple but 
lengthy. 

□ 

Example 2.3 If frequencies are Poisson distributed, ~ Poisson{X), then 

E[A^] = Var[A^] = E[(A^ - E[N])^] = A, 
E[(iV-E[iV])^] = A(1 + 3A), 



and compound loss moments calculated using Proposition 2.2 



are 



E[Z] = AE[Xi], Var[Z] = \E[Xf], E[{Z - E[Z]f] 
E[{Z - E[Z])^] = XE[Xt] + 3A2(E[X2])2. 



XE[X!], 



(19) 



Moreover, if the severities are lognormally distributed, Xi ~ CAf{fi, cr), then 

E[X'^] = exp{kn + Pay 2). (20) 



It is illustrative to see that in the case of compound Poisson, the moments 
can easily be derived using the following proposition. 

Proposition 2.4 (Cumulants of compound Poisson) The cumulants of 
the compound random variable Z = Xi + ■ ■ ■ + X^, where Xi, . . . , X^ are 
independent and identically distributed, and independent of N, are given by 

Kk = XE[X^], fc = l,2,... 



Proof 2 Using the definition of cumulants (17) and the characteristic func- 
tion for compound Poisson (Isj), calculate 



Inx(t) 
' dt^ 



k d'vit) 



t=o 



dt^ 



AE[Xf], A; = 1,2, 



i=0 



□ 



3 Value-at-Risk and Expected Shortfall 

Having calculated the compound loss distribution, the risk measures such as 
Value-at-Risk (VaR) and expected shortfall should be evaluated. Analyti- 
cally, VaR of the compound loss is calculated as the inverse of the compound 
distribution 

YaRa[Z] = H-\a) = mi{z e R : Pi[Z > z] < 1 - a} (21) 
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and the expected shortfall of the compound loss above the quantile qa 
VaRafZ], assuming that go > 0, is 



H{qo 



= E[Z\Z>qa] = 

nz] 

1 - H{qa) 1 - H{qa 



zh{z)dz 



1 



la 

qa 



zh{z)dz, 



(22) 



where E[Z] = E[A^]E[Xi] is the mean of compound loss Z. Note that ESq,[Z] 
is defined for a given quantile qa, that is, the quantile H~^{a) has to be com- 
puted first. It is easy to show (see formulas (40-43) in Luo and Shevchenko 
(2009)) that in the case of nonnegative severities, the above integral can be 
calculated via characteristic function as 



ESa[Z] 



X 



- Hiqa) 

E[Z] - H{qa)qo 



2qa 

IT 



Re [xix/qa)] 



cosx 



dx 



.(23) 



Remarks 3.1 



Strictly speaking, in the above formulas (22) and (23), we assumed that 
the quantile is positive, qa > 0, i.e. a > Pt[Z = 0] and we do not have 
complications due to discontinuity at zero. The case of go, = is not 
really important to operational risk practice, but can easily be treated 
if required. 



In the above formulas (22) and (23), H{qa) can be replaced by a. We 



kept H{qa), so that the formulas can easily be modified if expected 
exceedance E[Z\Z > L] should be calculated. In this case, should 
be replaced by L in these formulas. 

4 Monte Carlo Method 

The easiest numerical method to calculate the compound loss distribution is 
Monte Carlo (MC) with the following logical steps. 

Algorithm 4.1 (Monte Carlo for compound loss distribution) 

1. Fork = 1,...,K 



(a) Simulate the number of events N from the frequency distribution; 

(b) Simulate independent severities Xi, . . . , Xjq from the severity dis- 
tribution: 
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(c) Calculate Zk — Y^^^^Xi. 

2. Next k (i.e. do an increment k — k + 1 and return to step 1). 
All random numbers simulated in the above are independent. 

Obtained are samples from a compound distribution H{-). 

Distribution characteristics can be estimated using the simulated samples in 
the usual way described in many textbooks. Here, we just mention the quan- 
tile and expected shortfall which are of primary importance for operational 
risk. 



4.1 Quantile Estimate 

Denote samples Zi, . . . , Zk sorted into the ascending order as Zi < . . . < Zk, 
then a standard estimator of the quantile — H~^{a) is 

Qa = ZyKa\+l- (24) 

Here, [.J denotes rounding downward. Then, for a given realisation of the 

sample Z = the quantile estimate is = ZYKa\+i- It is important to 
estimate numerical error (due to the finite number of simulations K) in the 
quantile estimator. Formally, it can be assessed using the following asymp- 
totic result 

^^^"^'^:(g,-5,)^A^(0,l), as X^oo; (25) 



a] 



see e.g. Stuart and Ord (1994, pp.356-358) and Glasserman (2004, p.490). 
This means that the quantile estimator converges to the true value g„ as 
the sample size K increases and asymptotically Qa is normally distributed 
with the mean and standard deviation 

stdev[Q«] = W > . (26) 

However, the density h{qa) is not known and the use of the above formula 
is difficult. In practice, the error of the quantile estimator is calculated 
using a non-parametric statistic by forming a conservative confidence interval 
[Z'^'^\ to contain the true quantile value qa with the probability at least 

Pr[Zr <qa<Zs]>J, l<r <s<K. (27) 

Indices r and s can be found by utilising the fact that the true quantile qa 
is located between Zm and Zm+i for some M. The number of losses M not 
exceeding the quantile qa has a binomial distribution, Bin{K, a), because it 
is the number of successes from K independent and identical attempts with 
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success probability a. Thus the probabihty that the interval [Z^, contains 
the true quantile is simply 

Pr[r<M<s-l] = ^(' (28) 

i=r ^ ^ 

One typically tries to choose r and s that are symmetric around and closest 



to the index \Ka\ + 1, and such that the probability (28) is not less than 
the desired confidence level 7. The mean and variance of the binomial dis- 
tribution are Ka and Ka{\ — a) respectively. For large approximating 
the binomial by the normal distribution with these mean and variance leads 
to a simple approximation for the conservative confidence interval bounds: 



r = [/J, l = Ka- F^\{1 + ^)/2)^Ka{l - «), 

s = \u], u = Ka + F^\{1 + ^)/2)^Ka{l - a), (29) 

where [.] denotes rounding upwards and F^^{-) is the inverse of the standard 
normal distribution A/'(0, 1). The above formula works very well for Ka{l — 
a) > 50 approximately. 

Remarks 4.2 

• A large number of simulations, typically K > 10^, should be used to 
achieve a good numerical accuracy for the 0.999 quantile. However, a 
priori, the number of simulations required to achieve a specific accuracy 
is not known. One of the approaches is to continue simulations until a 
desired numerical accuracy is achieved. 

• If the number of simulations to get acceptable accuracy is very large 
(e.g. K > 10^) then you might not be able to store the whole array 
of samples Zi, . . . , Zk when implementing the algorithm, due to com- 
puter memory limitations. However, if you need to calculate just the 
high quantiles then you need to save only \Ka\ + 1 largest samples to 



estimate the quantile (24). This can be done by using the sorting on the 
fly algorithms, where you keep a specified number of largest samples 
as you generate the new samples; see Press et al (2002, Section 8.5). 
Moments (mean, variance, etc) can also be easily calculated on the fly 
without saving all samples into the computer memory. 



To use (29) for estimation of the quantile numerical error, it is im- 



portant that MC samples Zi, . . . , Zk are independent and identically 



distributed. If the samples are correlated, then (29) can significantly 
underestimate the error. In this case, one can use batch sampling or 
effective sample size methods; see e.g. Kass et al (1998). 

Example 4.3 Assume that if = 5 x 10^ independent samples were drawn 
from £A/'(0,2). Suppose that we would like to construct a conservative 
confidence interval to contain the 0.999 quantile with probability at least 



7 = 0.95. Then, sort the samples in ascending order and using (29) calculate 



F-i((l + 7)/2) ^ 1.96, r = 49936 and s = 49964 and [Ka\ + 1 = 49951 
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4.2 Expected Shortfall Estimate 



Given independent samples Zi, . . . , Zx from the same distribution and the 
estimator of VaRo[2'], a typical estimator for expected shortfall Ua = 
E[Z\Z > VaR„[Z]] is 



"~ 1 ^ ~ K-\Ka\ ■ ^ ' 



Here, 1{.} is a standard indicator symbol defined as 1 if condition in {•} 



is true and otherwise. Formula (30) gives an expected shortfall estimate 



Qa for a given sample realisation, Z = z. From the strong law of large 
numbers applied to the numerator and denominator and the convergence of 



the quantile estimator (25), it is clear that 

Via — )■ (jJa (31) 

with probability 1, as the sample size increases. If we assume that the quan- 
tile Qa is known, then in the limit K — )■ oo, the central limit theorem gives 



^Ar(0,l), (32) 



where a, for a given realisation Z = z, can be estimated as 

^2 _ r^Sfc=l(^fc ~ ^aY'^Zk>qa 



a' = K- 



1 ^Zk>qa 



Then, the standard deviation of Vt^ is estimated by ojyK; see Glasser- 
man (2005). However, it will underestimate the error in expected shortfall 
estimate because the quantile is not known and estimated itself by g^. Ap- 
proximation for asymptotic standard deviation of expected shortfall estimate 
can be found in Yamai and Yoshiba (2002, Appendix 1). In general, the stan- 
dard deviation of the MC estimates can always be evaluated by simulating 
K samples many times. For heavy-tailed distributions and high quantiles, it 
is typically observed that the error in quantile estimate is much smaller than 
the error in expected shortfall estimate. 

Remarks 4.4 Expected shortfall does not exist for distributions with infi- 
nite mean. Such distributions were reported in the analysis of operational 
risk losses; see Moscadelli (2004). 



5 Panjer Recursion 

It appears that, for some class of frequency distributions, the compound 
distribution calculation via the convolution (|4]) can be reduced to a simple 
recursion introduced by Panjer (1981) and referred to as Panjer recursion. A 
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good introduction of this method in the context of operational risk can be 
found in Panjer (2006, Sections 5 and 6). Also, a detailed treatment of Panjer 
recursion and its extensions is given in a recently published book Sundt and 
Vernic (2009). Below we summarise the method and discuss implementation 
issues. 

Firstly, Panjer recursion is designed for discrete severities. Thus, to apply 
the method for operational risk, where severities are typically continuous, the 
continuous severity should be replaced with the discrete one. For example, 
one can round all amounts to the nearest multiple of monetary unit 6, e.g. 
to the nearest USD 1000. Define 

/, = Pr[Xi = M], p, = FT[N = k], hk = FT[Z = k6], (33) 

with /o = and k = 0,1, . . . . Then, the discrete version of ^ is 



k=l 

ho = Pt[Z = 0]= Pt[N = 0]= Po, (34) 

where fi'> = Eto ft"^* with fl^> = 1 and fi'^* = if n > 1. 
Remarks 5.1 

• Note that the condition /o = Pr[Xi = 0] = implies that /n^''* = for 
k > n and thus the above summation is up to n only. 

• If /o > 0, then > for all n and k; and the upper limit in 
summation (34) should be replaced by infinity. 



• The number of operations to calculate ho, hi, . . . , hn using (34) explic- 
itly is of the order of n^. 

If the maximum value for which the compound distribution should be 
calculated is large, the number of computations become prohibitive due to 
0{n^) operations. Fortunately, if the frequency N belongs to the so-called 
Panjer classes, (34) is reduced to a simple recursion introduced by Panjer 
(1981) and referred to as Panjer recursion. 

Theorem 5.2 (Panjer recursion) If the frequency probability mass func- 
tion Pn, n = 0,1, . . . satisfies 

Pn = ( 0. -\ — ) Pn-i, for n > 1 and a, 6 G M, (35) 



then it is said to be in Panjer class (a, b, 0) and the compound distribution 
(34) satisfies the recursion 

1 f , bj 



hn = z ^y](«H ) fjhn-j, n>l, 

oo 

ho = 5^(/o)Vfe- (36) 



fe=0 
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The initial condition in (36) is simply a probability generating function 
of N at /o, i.e. ho = -ipifo), see If /o = 0, then it simplifies to ho = po- It 
was shown in Sundt and Jewell (1981), that (35) is satisfied for the Poisson, 
negative binomial and binomial distributions. The parameters (a, b) and 
starting values ho are listed in Table [T} 

Remarks 5.3 

• If severity is restricted by a value of the largest possible loss m, then 
the upper limit in the recursion (36) should be replaced by min(m,'n.). 



The Panjer recursion requires O(n^) operations to calculate ho, . . . ,hn 
in comparison with asymptotic 0{n^) of explicit convolution. 

Strong stability of Panjer recursion was established for the Poisson and 
negative binomial cases; see Panjer and Wang (1993). The accumulated 
rounding error of the recursion increases linearly in n with a slope not 
exceeding one. Serious numerical problems may occur for the case of 
binomial distribution. Typically, instabilities in the recursion appear 
for significantly underdispersed frequencies of severities with a large 
negative skewness which are not typical in operational risk. 

In the case of severities from a phase-type distribution (distribution 



with a rational probability generating function), the recursion (36) is 
reduced to 0(n) operations; see Hipp (2003). Typically, the sever- 
ity distributions are not phase-type distributions and approximation 
is required. This is useful for modelling small losses but not suitable 
for heavy-tailed distributions because the phase-type distributions are 
light tailed; see Bladt (2005) for a review. 

The Panjer recursion can be implemented as follows: 
Algorithm 5.4 (Panjer recursion) 

1. Initialization: calculate fo and ho, see Table\I^ and set Ho = ho- 

2. For n = 1,2,... 

(a) Calculate fn- If severity distribution is continuous, then fn can be 



found as described in Section \57l\ 

'"n—j ? 



(b) Calculate h^ = Yl]=i + ^) fj^^r 

(c) Calculate = -f^n-i + ^n/ 

(d) Interrupt the procedure if Hn is larger than the required quantile 
level a, e.g. a = 0.999. Then the estimate of the quantile g„ is 
n X 6. 

3. Next n (i.e. do an increment n = n + 1 and return to step 2). 
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5.1 Discretisation 



Typically, severity distributions are continuous and thus discretisation is re- 
quired. To concentrate severity, whose continuous distribution is F[x), on 
{0, 6, 26, . . .}, one can choose 6 > and use the central difference approxima- 
tion 

/o = F{S/2), 

fn = F{n6 + 6/2) - F{n6 - 6/2), n = 1, 2, . . . . (37) 

Then the compound discrete density /i„ is calculated using Panjer recursion 
and compound distribution is calculated as Hn = Y^^=o^i- an example, 
Table[2]gives results of calculation of the Pozssor;,(100) — £A/'(0, 2) compound 
distribution up to the 0.999 quantile in the case of step 6 = USD 1. Of course 
the accuracy of the result depends on the step size as shown by the results for 
the 0.999 quantile vs 6, see Table |3] and Figure [1} It is, however, important 
to note that the error of the result is due to discretisation only and there is 
no truncation error (i.e. the severity is not truncated by some large value). 
Discretisation can also be done via the forward and backward differences: 

= Fin6 + 6)-Fin6); fi^, = F{n6) - F{n6 ~ 6). (38) 

These allow for calculation of the upper and lower bounds for the compound 
distribution: 

n n 

K = Y.^'^: Hi: = Y,hi (39) 

1=0 1=0 

For example, see Table |4] presenting results for Poisson(lOO) — £A/'(0, 2) 
compound distribution calculated using central, forward and backward dif- 
ferences with step 6 = USDl. The use of the forward difference gives 
the upper bound for the compound distribution and the use of gives the 
lower bound. Thus the lower and upper bounds for a quantile are obtained 
with and respectively. In the case of Table |4] example, the quantile 
bound interval is [USD 5811, USD 5914] with the estimate from the central 
difference USD 5849. 



5.2 Computational Issues 

UnderfloAivj^ in computations of (36) will occur for large frequencies during 
the initialization of the recursion. This can easily be seen for the case of 
Poisson{X) and /o = when ho = exp(— A), that is, the underflow will 
occur for A > 700 on a 32bit computer with double precision calculations. 
Re-scahng ho by large factor 7 to calculate the recursion (and de-scahng the 

^Underflow/overflow are the cases when the computer calculations produce a number 
outside the range of representable numbers leading or ±00 outputs respectively. 
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result) does not really help because overflow will occur for '-yh{n). The follow- 
ing identity helps to overcome this problem in the case of Poisson frequency: 

ij(™)*(z;A/m) = /7(z;A). (40) 

That is, calculate the compound distribution H{z\ A/m) for some large m to 
avoid underflow. Then preform m convolutions for the obtained distribution 
directly or via FFT; see Panjer and Willmot (1986). Similar identity is 
available for negative binomial, NegBin{r,p): 

H^'^>{z;r/m) = H{z]r). (41) 
In the case of binomial, Bin{M,p): 

H^'"'>{z; mi) * H{z; m^) = H{z; M), (42) 

where mi = \_M/m\ and m2 = M — mim. 

For efficiency, one can choose m = 2^^ so that instead of m convolutions 
of H{-) only k convolutions are required if^^^*, H^'^^*, . . . , H^^ )*, where each 
term is the convolution of the previous one with itself. 



5.3 Panjer Extensions 



The Panjer recursion formula (36) can be extended to a class of frequency 
distributions (a, 6,1). 

Definition 5.5 (Panjer class (a, b, 1)) The distribution is said to be in 
(a, b, 1) Panjer class if it satisfies 

Pn = — j Pn-i, for n >2 and a, 6 G M. (43) 



Theorem 5.6 (Extended Panjer recursion) For the frequency distribu- 
tions in a class (a, b, 1): 

{pi-{a + b)po)fn + Y^'^=i{a + bj/n)fjhn-j 

hn = z -? , > 1, 

1 - afo 

oo 

ho = (44) 

fc=0 

The distributions of (a, 6, 0) class are special cases of (a, 6, 1) class. There 
are two types of frequency distributions in (a, 6, 1) class: 

• zero-truncated distributions, where po = 0: i.e. zero truncated Poisson, 
zero truncated binomial and zero-truncated negative binomial. 

• zero-modified distributions, where po > 0: the distributions of (a, 6, 0) 
with modified probability of zero. It can be viewed as a mixture of 
(a, 6, 0) distribution and degenerate distribution concentrated at zero. 



15 



Finally, we would like to mention a generalization of Panjer recursion for 
the (a, b, I) class 

Pn= (^a + ^ pn-i, for n > / + 1. (45) 

For initial values Po = ■ ■ ■ = pi-i = 0, and in the case of /o = 0, it leads to 
the recursion 

n 

K = Piff* + 5^ (a + bj/n) f,K_„ n>l. 
i=i 

The distribution in this class is, for example, / — 1 truncated Poisson. For an 
overview of high order Panjer recursions, see Hess et al (2002). Other types 
of recursions 

k 

Pn = "^{o-j + bj/n)pn^i, n>l, (47) 
i=i 

are discussed in Sundt (1992). Application of the standard Panjer recursion 
in the case of the generalised frequency distributions such as the extended 
negative binomial, can lead to numerical instabilities. Generalization of the 
Panjer recursion that leads to numerically stable algorithms for these cases 
is presented in Gerhold et al (2009). Discussion on multivariate version 
of Panjer recursion can be found in Sundt (1999) and bivariate cases are 
discussed in Vernic (1999) and Hesselager (1996). 



5.4 Panjer Recursion for Continuous Severity 

The Panjer recursion is developed for the case of discrete severities. The 
analog of Panjer recursion for the case of continuous severities is given by 
the following integral equation. 

Theorem 5.7 (Panjer recursion for continuous severities) For fre- 
quency distributions in (a, b, 1) class and continuous severity distributions on 
positive real line: 

h{z) = pifiz) + / (a + by/z)f{y)h{z - y)dy. (48) 
Jo 

The proof is presented in Panjer and Willmot (1992, Theorem 6.14.1 
and 6.16.1). Note that the above integral equation holds for (a, 6,0) class 
because it is a special case of (a, 6, 1). The integral equation ( [48] ) is a Volterra 
integral equation of the second type. There are different methods to solve it 
described in Panjer and Willmot (1992). A method of solving this equation 
using hybrid MCMC (minimum variance importance sampling via reversible 
jump MCMC) is presented in Peters et al (2007) . 
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6 Fast Fourier Transform 



The FFT is another efficient method to calculate compound distributions via 
the inversion of the characteristic function. The method has been known for 
many decades and originates from the signal processing field. The existence of 
the algorithm became generally known in the mid-1960s, but it was indepen- 
dently discovered by many researchers much earlier. One of the early books 
on FFT is Brigham (1974). A detailed explanation of the method in appli- 
cation to aggregate loss distribution can be found in Robertson (1992). In 
our experience, operational risk practitioners in banking regard the method 
as difficult and rarely use it in practice. In fact, it is a very simple algorithm 
to implement, although to make it really efficient, especially for heavy-tailed 
distribution, some improvements are required. Below we describe the essen- 
tial steps and theory required for successful implementation of the FFT for 
operational risk. 

As with Panjer recursion case, FFT works with discrete severity and based 
on the discrete Fourier transformation defined as follows. 

Definition 6.1 (Discrete Fourier transformation) 

For a sequence /o, /i, . . . , fu-i, the discrete Fourier transformation (DFT) 
is defined as 

M-l . . V 

^/mexpl ^mfcj , A; = 0,1,...,M-1 (49) 

m=0 ^ ^ 

and the original sequence fk can be recovered from (pk by the inverse trans- 
formation 

1 ^^'^ ( 2'Ki \ 

fk=MJ2'f'rnexp(- -^mk j , A; = 0, 1, . . . , M - 1. (50) 



Here, M is some truncation point. It is easy to see that to calculate M 
points of (f)m: the number of operations is of the order of M^, i.e. O(M^). If 
M is a power of 2, then DFT can be efficiently calculated via FFT algorithms 
with the number of computations 0(Af log2 M). This is due to the property 
that DFT of length M can be represented as the sum of DFT over even 
points (pl and DFT over odd points 0^: 

/ 2TTi \ 

(f)k = (f)l + expl—k](Pl] 

M/2-1 

^ Yl /2m exp 
m=0 
M/2-1 

<^fe = X] /2m+iexp 

m=0 
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Subsequently, each of these two DFTs can be calculated as a sum of two 
DFTs of length M/4. For example, (pl is calculated as a sum of (pl^ and 0^°. 
This procedure is continued until the transforms of the length 1. The latter 
is simply identity operation. Thus every obtained pattern of odd and even 
DFTs will be fm for some m: 

ieO'"Ooe r 

(Pk — Jm- 

The bit reversal procedure can be used to find m that corresponds to a 
specific pattern. That is, set e = and o = 1, then the reverse pattern of 
e's and o's is the value of m in binary. Thus the logical steps of FFT are as 
follows. 

Algorithm 6.2 (Simple FFT) 

1. Sort the data in a hit-reversed order. The obtained points are simply 
one-point transforms. 

2. Combine the neighbor points into non- overlapping pairs to get two-point 
transforms. Then combine two-point transforms into 4-point trans- 
forms and continue subsequently until the final M point transform is 
obtained. Thus there are log2 M iterations and each iteration involves 
of the order of M operations. 

The implementation of a basic FFT algorithm is very simple; correspond- 
ing C or Fortran codes can be found in Press et al (2002, Chapter 12). 

6.1 Compound Distribution via FFT 

Calculation of the compound distribution via FFT can be done using the 
following logical steps. 

Algorithm 6.3 (Compound Distribution via FFT) 

1. Discretise severity to obtain 

where M = T' with integer r and M is the truncation point in the 
aggregate distribution; 

2. Using FFT, calculate the characteristic function of the severity 

'^0, ■ ■ ■ , V^M-i; 

3. Calculate the characteristic function of the compound distribution using 

Xm = ipi'fm), m = 0, 1, . . . ,M - 1. 
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4- Perform inverse FFT (which is the same as FFT except the change of 
sign under the exponent and factor l/M) applied to xo, • • • ,Xm-i to 
obtain the compound distribution ho,hi, . . . , Hm-i- 

Remarks 6.4 To calculate the compound distribution in the case of the 
severity distribution F{x) with a finite support (i.e. 0<a<a:<6< oo) 
one can set F{x) = for x outside the support range when calculating 



discretised severity /o,...,/m-i using (37). For example, this is the case 
for distribution of losses exceeding some threshold. Note that we need to 
set F[x) = in the range x G [0, a) due to the finite probability of zero 
compound loss. 



6.2 Aliasing Error and Tilting 

If there is no truncation error in the severity discretisation, i.e. Ylm=o fm = ^, 
then FFT procedure calculates the compound distribution on m = 0, 1, . . . , M. 
That is, the mass of compound distribution beyond M is "wrapped" and ap- 
pears in the range m = 0, . . . , M — 1 (the so-called aliasing error). This error 
is larger for heavy-tailed severities. To decrease the error for compound dis- 
tribution on 0, 1, . . . , n, one has to take M much larger than n. If the severity 
distribution is bounded and M is larger than the bound, then one can put 
zero values for points above the bound (the so-called padding by zeros). An- 
other way to reduce the error is to apply some transformation to increase 
the tail decay (the so-called tilting). The exponential tilting technique for 
reducing aliasing error under the context of calculating compound distribu- 
tion was first investigated by Grubel and Hermesmeier (1999). Many authors 
suggest the following tilting transformation: 

= exp(-j^)/„ j = 0,l,...,M-l, (51) 

where 6 > 0. This transformation commutes with convolution in a sense 
that convolution of two functions f{x) and g{x) equals the convolution of 
the transformed functions f{x) = f{x)exp{—6x) and g{x) = g{x) exp{—6x) 
multiplied by exp(6'x), i.e. 

{f*g){x) = e'%f*g){x)- (52) 

This can easily be shown using the definition of convolution. Then calculation 
of the compound distribution is performed using the transformed severity 
distribution as follows. 

Algorithm 6.5 (Compound distribution via FFT with tilting) 

1. Define /o, /i, . . . , /m-i for some large M; 

2. Perform tilting, i.e. calculate the transformed function fj = exp{—j6)fj, 
j = 0,l,...,M-l; 

3. Apply FFT to a set /o, . . . , fhi-i to obtain (pQ, . . . , 0a/-i; 
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4- Calculate Xm = 4'{4>m),m = 0, 1, . . . , M — 1; 

5. Apply the inverse FFT to the set xo, • • • ? Xm-i, to obtain Hq, . . . , Hm-i; 

6. Untilt by calculating final compound distribution as hj = hj exp{6j) . 

This tilting procedure is very effective in reducing the ahasing error. The 
parameter 9 should be as large as possible but not producing under- or over- 
ffow that will occur for very large 6. It was reported in Embrechts and Frei 
(2009) that the choice M6 ~ 20 works well for standard double precision (8 
bytes) calculations. Evaluation of the probability generating function ip{-) of 
the frequency distribution may lead to the problem of underflow in the case 
of large frequencies that can be resolved using methods described in Section 



Example 6.6 To demonstrate the effectiveness of the tilting, consider the 
following calculations: 

• FFT with the central difference discretisation, where the tail proba- 
bility compressed into the last point /a/-i = 1 — F{6{M — 1) — 6/2). 



(1) 



Denote the corresponding quantile estimator as Qo 999, 

• FFT with the central difference discretisation with the tail probability 
ignored, i.e. /m-i = F{6{M - 1) + 6/2) - F{6{M - 1) - 6/2). Denote 
the corresponding quantile estimator as Q0.999! 

• FFT with the central difference discretisation utilising tilting Q^ggg. 
The tilting parameter 6 is chosen to be 6* = 20/M. 

The calculation results presented in Table |5] demonstrate the efficiency 
of the tilting. If FFT is performed without tilting then the truncation level 
for the severity should exceed the quantile significantly. In this particular 
case it should exceed by approximately factor of 10 to get the exact result 
for this discretisation step. The latter is obtained by Panjer recursion that 
does not require the discretisation beyond the calculated quantile. Thus the 
FFT and Panjer recursion are approximately the same in terms of computing 
time required for quantile estimate in this case. However, once the tilting is 
utilised, the cut off level does not need to exceed the quantile significantly to 
obtain the exact result - making FFT superior to Panjer recursion. In this 
example, the computing tim^for FFT with tilting is 0.17sec in comparison 
with 5.76sec of Panjer recursion, see Table|3| Also, in this case, the treatment 
of the severity tail by ignoring it or absorbing into the last point fpi-i does 
not make any difference when tilting is applied. 



^Computing time is quoted for a standard Dell laptop Latitude D820 with Intel(R) 
CPU T2600 @ 2.16 GHz and 3.25 GB of RAM. 
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7 Direct Numerical Integration 



In the case of nonnegative severities, the distribution of the compound loss 
is given by (12), i.e. 

oo 

H{z) = - [ Re[x{t)f^^dt, z > 0, (53) 

IT J t 


where x(^) is a compound distribution characteristic function calculated via 
the severity characteristic function ip(t) using ([t]). For example, the explicit 
expression of Re[x(t)] for Poisson{X) is 

Re[x(t)] = e-^exp(ARe[¥?(t)]) x cos(AIm[(/?(t)]). (54) 
Hereafter, direct calculation of the distribution function for annual loss Z 



using (53) is referred to as direct numerical integration (DNI). 

Much work has been done in the last few decades in the general area 
of inverting characteristic functions numerically. Just to mention a few, see 
the works by Bohman (1975); Seal (1977); Abate and Whitt (1992), (1995); 
Heckman and Meyers (1983); Shephard (1991); Waller et al (1995); and Den 
Iseger (2006). These papers address various issues such as singularity at 
the origin; treatment of long tails in the infinite integration; and choices 
of quadrature rules covering different objectives with different distributions. 
Craddock et al (2000) gave an extensive survey of numerical techniques for 
inverting characteristic functions. 

Each of the many existing techniques has particular strengths and weak- 
nesses, and no method works equally well for all classes of problems. In 
an operational risk context, for instance, there is a special need in com- 
puting the 0.999 quantile of the aggregate loss distribution. The accuracy 
demanded is high and at the same time the numerical inversion could be very 
time consuming due to rapid oscillations and slow decay in the characteristic 
function. This is the case, for example, for heavy-tailed severities. Also, 
the characteristic function of compound distributions should be calculated 
numerically through semi-infinite integrations. A tailor-made numerical al- 



gorithm to integrate (53) was presented in Luo et al (2007) and Luo and 
Shevchenko (2009) with a specific requirement on accuracy and efficiency in 
calculating high quantiles such as 0.999 quantile. The method works well 
for both a wide range of frequencies from very low to very high (> 10^) and 
heavy-tailed severities. 



7.1 Forward and Inverse Integrations 

The task of the characteristic function inversion is analytically straightfor- 
ward, but numerically difficult in terms of achieving high accuracy and com- 
putational efficiency simultaneously. 

Accurate calculation of the high quantile as an inverse of the distribution 
function requires high precision in evaluation of the distribution function. 
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To demonstrate, consider the lognormal distribution £A/'(0,2). In this case, 
the "exact" 0.999 quantile go.999 = 483.2164. . . . However, at a = 0.99902, 
the quantile becomes = 489.045 .... That is, a mere 0.002% change in 
the distribution function value causes more than 1% change in the quantile 
value. In the case of a compound distribution, the requirement for accuracy 
in the distribution function could be even higher, because 1/ f{x) could be 
larger at x = go.999- Note that, the error propagation from the distribution 
function level to the quantile value is implied by the relation between the 
density /(x) and its distribution function F{x): dF/dx = f{x). 

The computation of compound distribution through the characteristic 
function involves two steps: computing the characteristic function (Fourier 
transform of the density function, referred to as the forward integration) and 
inverting it (referred to as the inverse integration). 



7.1.1 Forward Integration 

This step requires integration (|5|, that is, calculation of the real and imagi- 
nary parts of the characteristic function for a severity distribution: 

00 00 

Re[ip{t)] = J f{x) cos{tx)dx, lm[Lp{t)] = J f{x)sm{tx)dx. (55) 



Then, the characteristic function of the compound loss is calculated using 
([7]). These tasks are relatively simple because the severity density typically 
has closed-form expression, and is well-behaved having a single mode. 

This step can be done more or less routinely and many existing algorithms 
can be employed. The oscillatory nature of the integrand only comes from the 
sin( ) or cos( ) functions. This well-behaved weighted oscillatory integrand 
can be effectively dealt with by the modified Clenshaw- Curtis integration 
method; see Clenshaw and Curtis (1960) and Piessens et al (1983). In this 
method the oscillatory part of the integrand is transferred to a weight func- 
tion, the non-oscillatory part is replaced by its expansion in terms of a finite 
number of Chebyshev polynomials and the modified Chebyshev moments are 
calculated. If the oscillation is slow when the argument t of the characteris- 
tic function is small, the standard Guass-Legendre and Kronrod quadrature 
formulae are more effective; see Kronrod (1965), Golub and Welsh (1969), 



Szego (1975), and Section 7.2 In general, double precision accuracy can 
be routinely achieved for the forward integrations using standard adaptive 
integration functions commonly available in many software packages. 



7.1.2 Inverse Integration 

This step requires integration (53), which is much more challenging task 
Changing variable x = t x z, (53) can be rewritten as 

2 Re[x{x/z)] 

TT X ' 



H(z) 



G{x, z) sm{x)dx, G{x,z) 



(56) 
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where x{t) depends on Re[v9(t)] and Im[y9(t)] calculated from the forward 
semi-infinite integrations ( [55] ) for any required argument t. The total num- 
ber of forward integrations required by the inversion is usually quite large. 
This is because in this case the characteristic function could be highly os- 
cillatory due to high frequency and it may decay very slowly due to heavy 
tails. There are two oscillatory components in the integrand represented by 
sin(x) and another part in Re[x(a;/-2)]. It is convenient to treat sin(x) as 
the principal oscillatory factor and the other part as secondary. Typically, 
given z, Re[x{x/z)] decays fast initially and then approaches zero slowly as 
X approaches infinity. 



To calculate (56), one could apply the same standard general purpose 
adaptive integration routines as for the forward integration. However, this is 
typically not efficient because it does not address irregular oscillation specifi- 
cally and can lead to an excessive number of integrand evaluations. A simple 



approach that can be taken is to divide the integration range of (56) into 



intervals of equal length vr (referred to as vr-cycle) and truncate at 2K7r: 

H{z) ^ J2 ^k, Hk= / G{x) sm{x)dx. (57) 

k=0 L 



Within each vr-cycle, the secondary oscillation could be dominating for some 
early cycles, thus the vr-cycle could in fact contain multiple cycles due to 
the "secondary" oscillation. Thus a further sub-division is warranted. Sub- 
dividing interval (fcvr, {k + l)vr) into segments of equal length of = 
vr/rifc, (57) can be written as 



where 



afcj = fcvr + (j - l)Afc, hk,j = ak,j + Afc. 

The above calculation will be most effective if the sub-division is made adap- 
tive for each vr-cycle according to the changing behaviour of G{x). Assum- 
ing that for the first vr-cycle {k = 0) we have initial partition no, Luo and 
Shevchenko (2009) recommends making Uk adaptive for the subsequent cycles 
by the following two simple rules: 

• Let Hk be proportional to the number of vr-cycles of the secondary 
oscillation - the number of oscillations in G{x) within each principal 
vr-cycle; 

• Let Hk be proportional to the magnitude of the maximum gradient of 
G{x) within each principal vr-cycle. 

Application of these rules requires correct counting of secondary cycles and 
good approximation of the local gradient in G{x). Both can be achieved with 
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a significant number of points at wliicli G{x) is computed witliin eacli cycle 
using, for example, the m-point Gaussian quadrature described in tlie next 
section. 



7.2 Gaussian Quadrature for Subdivisions 

Witli a proper sub-division, even a simple trapezoidal rule can be applied to 
get a good approximation for integration over the sub-division Hjf in ( [58| . 
However, higher order numerical quadrature can achieve higher accuracy for 
the same computing effort or it requires less computing effort for the same 
accuracy. The m-point Gaussian quadrature makes the computed integral 
exact for all polynomials of degree 2m— 1 or less. In particular: 

/ g{x)dx^ -Y^w,g{{a + h + CA)/2), (59) 
^ i=i 

where < < 1 and — 1 < < 1 are the i^^ weight and the i*^ abscissa of 
the Gaussian quadrature respectively, A = 6 — a and m is the order of the 
Gaussian quadrature. 

Typically, even a simple 7-point Gaussian quadrature (m = 7), which 
calculates all polynomials of degree 13 or less exactly, can successfully be 



used to calculate ifp'* in (|57l 581). For completeness. Table |6| presents 7-point 



Gaussian quadrature weights and abscissas; other quadratures can be found 
in Piessens et al (1983). 

The efficiency of the Gaussian quadrature is much superior to the trape- 
zoidal rule. For instance, integrating the function sin(3a;) over the interval 
(0,7r), the 7-point Gaussian quadrature has a relative error less than 10~^, 
while the trapezoidal rule requires about 900 function evaluations (grid spac- 
ing 5x = 7r/900) to achieve a similar accuracy. The reduction of the number 



of integrand function evaluations is important for a fast integration of (57), 
because the integrand itself is a time consuming semi-infinite numerical in- 
tegration. 

The error of the m-point Gaussian quadrature rule can be accurately esti- 
mated if the 2m order derivative of the integrand can be computed (Kahaner 
et al (1989); Stoer and Bulirsch (2002)). In general, it is difficult to estimate 
the 2m order derivative and the actual error may be much less than a bound 
established by the derivative. As it has already been mentioned, a common 
practice is to use two numerical evaluations with the grid sizes different by 
the factor of two and estimate the error as the difference between the two 
results. Equivalently, different orders of quadrature can be used to estimate 
error. Often, Guass-Kronrod quadrature is used for this purpose. Adaptive 
integration functions in many numerical software packages use this estimate 
to achieve an overall error bound below the user-specified tolerance. 
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7.3 Tail Integration 



The truncation error of using (57) is 



oo 

Ht= j G{x) sm{x)dx. (60) 

2KiT 

For higher accuracy, instead of increasing truncation length at the cost of 
computing time, one can try to calculate the tail integration approxi- 



mately or use tilting transform (51). Integration of (60) by parts gives 

°? fe-i 

/ G{x)sm{x)dx = G(2ir7r) + ^(-l)JG(2^)(2ir7r) 



2Ktt 

"OO 

\k 



G^^''\x)sm{x)dx, (61) 



2K-K 



where k > I, G^'^^\2K'k) is the 2j-th order derivative of G{x) at the trunca- 
tion point. Under some conditions, as if — )■ oo, 

/ G{x) sm{x)dx G{2Kn) + ^(-l)^G(2i)^2if7r). 



2Kn 



For example, if we assume that for some 7 < 0, G^'^\x) = 0(x'^~"^), m = 
0, 1, 2, . . . as if — )■ 00, then the series converges to the integral. However, 
this is not true for some functions, such as exp(— x); typically in this case the 
truncation error is not material. It appears that often, the very first term in 
(61) gives a very good approximation 

00 

Ht= j G{x) sin(x)rfx ^ G{2Kn) (62) 

2Ktt 

for the tail integration or does not have a material impact on the overall 
integration; see Luo and Shevchenko (2009, 2010). This elegant result means 
that we only need to evaluate the integrand at one single point x = 2ttK for 



the entire tail integration. Thus the total integral approximation (57) can 
be improved by including tail correction giving 



2K~l 



H{z) ^ ^ i^fc + G(2iV7r). (63) 



k=0 



Remarks 7.1 The approximation (62) can be improved by including further 



terms if derivatives are easy to calculate, e.g. Ht ~ G{2Kti) — G^'^\2Kti). 
If the oscillating factor is cos(a;) instead of sin(a;), one can still derive a one- 



point formula similar to (61 ) by starting the tail integration at (2K — l/2)7r 
instead of 2Kit. 
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Of course there are more elaborate methods to treat the truncation error 
which are superior to a simple approximation (62) in terms of better accu- 
racy and broader applicability, such as some of the extrapolation methods 
proposed in Wynn (1956), Sidi (1980) and Sidi (1988). 



7.4 Error Sources and Numerical Example 

Table [T] shows the convergence of DNI results (seven digits), for truncation 
lengths 2 < i^' < 80 in the cases of tail correction included and ignored. 
One can see a material improvement from the tail correction. Also, as the 
truncation length increases, both estimators with the tail correction and 
without converge. In this particular case we calculate compound distribution 
Poisson{100)-CJ\f{0,2) at the level z = 5853.1. The latter is the value 
that corresponds to the 0.999 quantile (within 1st decimal place) of this 
distribution as has already been calculated by Panjer recursion; see Table |3] 
Of course, to calculate the quantile at the 0.999 level using DNI, a search 
algorithm such as bisection should be used that will require evaluation of 
distribution function many times (of the order of 10) increasing computing 
time. Comparing this with Tables |3] and [5} one can see that for this case 
DNI is faster than Panjer recursion while slower than FFT (with tilting) by 
a factor of 10. 

The final result of the inverse integration has three error sources: the 
discretisation error of the Gauss quadrature; the error from the tail approx- 
imation; and the error propagated from the error of the forward integration. 
These were analysed in Luo and Shevchenko (2009). It was shown that the 
propagation error is proportional to the forward integration error bound. At 
the extreme case of A = 10^, a single precision can still be readily achieved if 
the forward integration has a double precision. For very large A, the propa- 
gation error is likely the largest among the three error sources. Though some 
analytic formulas for error bounds are available, these are not very useful 
in practise because high order derivatives are involved, which is typical for 
analytical error bounds. An established and satisfactory practice is to use 
finer grids to estimate the error of the coarse grids. 



8 Comparison of Numerical Methods 

For comparison purposes, Tables[8]and|9]present results for the 0.999 quantile 
of compound distributions Poisson{X)-CAf{0,2) and Poisson{\)-GPD{l, 1) 
(with A = 0.1, 10, 10^), calculated by the DNI, FFT, Panjer and MC methods. 
Note that, with the shape parameter ^ = 1, GPD{^, f3) has infinite mean 
and all higher moments. For DNI, FFT and Panjer recursion methods, the 
results, accurate up to 5 significant digits, were obtained as follows: 

• For DNI algorithm we start with a relatively coarse grid {hq = 1) 
and short truncation length K = 25, and keep halving the grid size 
and doubling the truncation length until the difference in the 0.999 
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quantile is within required accuracy. The DNI algorithm computes 



distribution function, H{z), for any given level z by (53), one point at 



a time. Thus with DNI we have to resort to an iterative procedure to 



inverse (53). This requires evaluating (53) many times depending on 
the search algorithm employed and the initial guess. Here, a standard 
bisection algorithm is employed. Other methods (MC, Panjer recursion 
and FFT) have the advantage that they obtain the whole distribution 
in a single run. 

• For Panjer recursion, starting with a large step (e.g. 6 = 8) the step 
5 is successively reduced until the change in the result is smaller than 
the required accuracy. 

• For FFT with tilting, the same step 6 is used as the one in the Panjer 
recursion. If we would not know the Panjer recursion results, then 
we would successively reduce the step S (starting with some large step) 
until the change in the result is smaller than the required accuracy. The 
truncation length M = 2^ has to be large enough so that 6M > Qg is 
satisfied. We use the smallest possible integer r that allows to identify 
the quantile, typically such that 6M ^ 2Qg. Here, Qg is the quantile 
to be computed, which is not known a priori and some extra iteration 
is typically required. Also, the tilting parameter is set to = 20/M. 

• For the MC estimates, the number of simulations, Nmc (denoted by 
K in Section |4]), ranges from 10® to 10^, so that calculations are ac- 
complished within ^ 10 min. The error of the MC estimate is approx- 
imately proportional to 1/ a/ Nmc and the calculation time is approxi- 
mately proportional to Nmc- Thus the obtained results allow to judge 
how many simulations (time) is required to achieve a specific accuracy. 

The agreement between FFT, Panjer recursion and DNI estimates is per- 
fect. Also, the difference between these results and corresponding MC esti- 
mates is always within the two MC standard errors. However, the CPU time 
is very different across the methods: 



• The quoted CPU time for the MC results is of the order 10 min. How- 
ever, it is clear from the standard error results (recalling that the error 
is proportional to I/^/Nmc) that the CPU time, required to get the 
results accurate up to five significant digits, would be of the order of 
several days. Thus MC is the slowest method. 

• Typically, the CPU time for both Panjer recursion and FFT increase 
as A increases, while CPU time for DNI does not change significantly. 

• FFT is the fastest method, though at very high frequency A = 10"^, DNI 
performance is of a similar order. As reported in Luo and Shevchenko 
(2009), DNI becomes faster than FFT for higher frequencies A > 10^. 
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• Panjer recursion is always slower than FFT. It is faster than DNI for 
small frequencies and much slower for high frequencies. 

Finally note that, the FFT, Panjer recursion and DNI results were ob- 
tained by successive reduction of grid size (starting with a coarse grid) until 
the required accuracy is achieved. The quoted CPU time is for the last iter- 
ation in this procedure. Thus the results for CPU time should be treated as 
indicative only. For comparison of FFT and Panjer, also see Embrechts and 
Frei (2009), and Biihlmann (1984). 



9 Closed- Form Approximation 

There are several well-known approximations for the compound loss distribu- 
tion. These can be used with different success depending on the quantity to 
be calculated and distribution types. Even if the accuracy is not good, these 
approximations are certainly useful from the methodological point of view 
in helping to understand the model properties. Also, the quantile estimate 
derived from these approximations can successfully be used to set a cut-off 
level for FFT algorithms that will subsequently determine the quantile more 
precisely. 



9.1 Normal and Translated Gamma Approximations 

Many parametric distributions can be used as an approximation for a com- 
pound loss distribution by moment matching. This is because the moments 
of the compound loss can be calculated in closed-form. In particular, the 



first four moments are given in Proposition 2.2 Of course these can only 



be used if the required moments exist which is not the case for some heavy- 
tailed risks with infinite moments. Below we mention normal and translated 
gamma approximations, discussed e.g. in McNeil et al (2005, Section 10.2.3). 



9.1.1 Normal Approximation 

As the severities Xi,X2, . . . are independent and identically distributed, at 
very high frequencies the central limit theory is expected to provide a good 
approximation to the distribution of the annual loss Z (if the second moment 
of severities is finite). Then the compound distribution is approximated by 
the normal distribution with the mean and variance given in Proposition |2.2[ 
that is, 

E{z) ^ M{Y\Zl VVar[Z]). (64) 

This is an asymptotic result and a priori we do not know how well it will 
perform for a specific distribution types and distribution parameter values. 
Also, it cannot be used for the cases where variance or mean are infinite. 
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Example 9.1 If N is distributed from Poisson{X) and Xi, . . . are in- 
dependent random variables from £A/'(/i, cr), then 

E[Z] = Aexp(/i + 0.5(j2), Var[Z] = Aexp(2/i + 2a'^). (65) 
9.1.2 Translated Gamma Approximation 



From (19), the skewness of the compound distribution, in the case of Poisson 



distributed frequencies, is 

E[(z-E[zi)3i mx^\ , , 

(Var[Z])'/' (AE[X2])^/2 
that approaches zero as A increases but finite positive for finite A > 0. To 



improve the normal approximation (64), the compound loss can be approxi 



mated by the shifted gamma distribution which has a positive skewness, that 
is, Z is approximated as y + a where a is a shift and F is a random variable 
from Gamma{a, [}). The three parameters are estimated by matching the 
mean, variance and skewness of the approximate distribution and the correct 
one: 

a + a/3 = E[Z]; a/?^ = Var[Z]; 4= = - E[Z])3]/ (Var[Z])^/^ (67) 



a 

This approximation requires the existence of the first three moments and 
thus cannot be used if the third moment does not exist. 

Example 9.2 If frequencies are Poisson distributed, ~ Poisson{X), then 
a + a/3 = XE[X]; a^^ = XE\X\ = XE\X^]/ (XE\X^]f^\ (68) 



9.2 VaR Closed-Form Approximation 

If severities Xi, . . . ,X^ are independent and identically distributed from the 
sub-exponential (heavy tail) distribution F{x), and frequency distribution 
satisfies 

oo 

Y^{l + eYVi[N = n] < oo 

n=0 

for some e > 0, then the tail of the compound distribution H{z), of the 
compound loss Z = Xi + ■ ■ ■ + Xjy, is related to the severity tail as 

1- H{z) ^E[N]{1- F{z)), as z ^ oo; (69) 

see Theorem 1.3.9 in Embrechts et al (1997). The validity of this asymptotic 
result was demonstrated for the cases when N is distributed from Poisson, 
binomial or negative binomial. It can be used to find the quantile of the 
annual loss 

VaR,[Z] ^ (^1 - , as « ^ 1. (70) 
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For application in the operational risk context, see Bocker and Kliippelberg 
(2005). Under the assumption that the severity has a finite mean, Bocker 
and SprittuUa (2006) derived a correction reducing the approximation error 



of (70) 



Example 9.3 Consider a heavy-tailed compound distribution Poisson{X)- 



GPD{i,(5). In this case, (70) gives 



VaR„[Z] ^ ^ ( ) ' as a ^ I. (71) 



This implies a simple scaling, VaRQ,[Z'] oc A^, with respect to the event in- 
tensity A for large a. 

Example 9.4 To demonstrate the accuracy the above approximations, con- 
sider compound distribution Poisson{X = 100)-CAf{fi = 0,a = 2) with rela- 
tively heavy tail severity. Calculating moments of the lognormal distribution 



E[X'"] using (20) and substituting into (19) gives 



E[Z] ^ 738.9056, Var[Z] ^ 298095.7987, 
E[(Z - E[Z])V(Var[Z])^/^ ^ 40.3428. 

Approximating the compound distribution by the normal distribution with 
these mean and variance gives normal approximation. Approximating the 



compound distribution by the translated gamma distribution (67) with these 
mean, variance and skewness gives: a ~ 0.002457, (3 ^ 11013.2329, a ^ 
711.8385. Figure |2^ shows the normal and translated gamma approxima- 
tions for the tail of the compound distribution. These are compared with 



the asymptotic result for heavy tail distributions (69) and "exact" values ob 



tained by FFT. It is easy to see that the heavy tail asymptotic approximation 



(69) converges to the "exact" result for large quantile level a — j- 1, while the 
normal and gamma approximations perform badly. The results for the case 
of not so heavy tail, when the severity distribution is CJ\f{0, 1), are shown in 
Figure [2]d. Here, the gamma approximation outperforms normal approxima- 
tion and heavy tail approximation is very bad. The accuracy of the heavy 



tail approximation (69) improves for more heavy-tailed distributions, such as 



GPD with infinite variance or even infinite mean. 



10 Conclusions 

In this paper we reviewed methods that can be used to calculate the dis- 
tribution of the compound loss. Overall, FFT with tilting is typically the 
fastest method though it involves tuning of the cut-off level, tilting param- 
eter and discretisation step. The easiest to implement is Panjer recursion 
that involves discretisation error only. DNI method is certainly competitive 
with FFT and Panjer for large frequencies, though its implementation can 
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be quite involved. Monte Carlo method is slow but simple in implementa- 
tion and it can easily handle multiple risks with dependence. The latter is 
problematic for FFT and Panjer recursion methods. In general, each of the 
reviewed techniques has particular strengths and weaknesses that a modeller 
should be aware of. The choice of the method is dictated by the specific 
objectives to be achieved. 

A List of Distributions 

Poisson distribution, Poisson{\). A Poisson distribution function is de- 
noted as Poisson{X). The random variable N has a Poisson distribution, 
denoted N ~ Poisson{X), if its probability mass function is 

\ k 

Pk = Pr[7V = A;] = -e-^ A > 0, A; e {0, 1, 2, . . .}. (72) 
Expectation, variance and variational coefficient are 

E[N] = A, Var[iV] = A, Vco[iV] = (73) 

V A 

Binomial distribution, Bin{n,p). A binomial distribution function is de- 
noted as Bin{n,p). The random variable has a binomial distribution, 
denoted N ~ Bin{n,p), if its probability mass function is 

p, = Fr[N = k] = (^l^ p\l - pr-\ p e (0,1), n e 1,2,. . . (74) 

for all k e {0, 1, . . . ,n}. Expectation, variance and variational coefficient are 



E[N] = np, VarlN] = np(l - p), Yco\N] = J- — -. (75) 

V np 

Negative binomial distribution, NegBin{r,p). A negative binomial dis- 
tribution function is denoted as NegBin{r,p). The random variable A^ has a 
negative binomial distribution, denoted A^ ~ NegBin{r,p), if its probability 
mass function is 

= Fr[N = A;] = ^^±^p^{l - p)\ p e (0, 1), r e (0, oo) (76) 

for all k G {0,1,2,...}. Here, V{r) is the gamma function. Expectation, 
variance and variational coefficient are 



E[Ar] = Var[Ar] = 'SL^l^ Yco[N] = (77) 

P P yjr{l-p) 
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Normal distribution, A/'(yU,o"). A normal (Gaussian) distribution function 
is denoted as J\f{fi,(r). The random variable X has a normal distribution, 
denoted X ~ a), if its probability density function is 

= Tib'"" (-^^) . > 0. " e « (78) 

for all X e M. Expectation, variance and variational coefficient are 

E[X] = /i, Var[X] = a^, Vco[X] = a/^. (79) 

Lognormal distribution, £N{ij,, a) A lognormal distribution function is 
denoted as jCN{iJ,,a). A random variable X has a lognormal distribution, 
denoted X ~ jCJ\f{ii, a), if its probability density function is 

for X > 0. Expectation, variance and variational coefficient are 



E[A:] = e'^+^"', Var[A:] = e2^+"'(e"' - 1), Vco[A:] = \/e-' - 1. (81) 

Gamma distribution, Gamma{a, /3). A gamma distribution function is 

denoted as Gamma{a, /3). The random variable X has a gamma distribution, 
denoted as X ~ Gamma{a, (3), if its probability density function is 

fix) = ^^^M-x/P), a > 0, /3 > (82) 
for X > 0. Expectation, variance and variational coefficient are 

E[X] = aP, Var[X] = Vco[X] = l/^/a. (83) 

Generalised Pareto distribution, GPD(^, (3). The GPD distribution 
function is denoted as GPD{C,, P). The random variable X has GPD distri- 
bution, denoted as X ~ GPD{^, P), if its distribution function is 

l-exp(-x//3), e = 0, ^^^^ 

where x > when ^ > and < x < —/3/^ when ^ < 0. The moments of 
X ~ GPD{^, /5), C ^ 0) can be calculated using 

/9"t7I 1 
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Table 1: Panjer recursion starting values ho and (a, b) parameters for Poisson, 
binomial and negative binomial distributions. 

a b ho 

Poisson{\) A exp(A(/o — 1)) 



NegBinir,q) 1-q (1 - q)(r - 1) (^l + (l _ 



Bin{rn, q) 


? 

1-q 


^ (i+g(/o-i)r 


Table 2: 


Example of Panjer 


recursion calculating the Poisson( 100) — 


£A/'(0, 2) compound distributions using central difference discretisation with 


the step 5 ■ 


- 1. 




n 


fn 


hn Hn 





0.364455845 


2.50419 X 10-28 2.50419 x lO-^^ 


1 


0.215872117 


5.40586 X 10-27 5.65628 x 10-^^ 


2 


0.096248034 


6.07589 X 10-26 6.64152 x lO-^^ 


5847 


2.81060 X 10-9 


4.44337 X 10-7 0.998999329 


5848 


2.80907 X 10-9 


4.44061 X 10-7 0.998999773 


5849 


2.80755 X 10-9 


4.43785 X 10-^ 0.999000217 
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Table 3: 


Convergence of Panjer recursion estimate, $0,999; of the 0.999 quan- 


tile for the Poisson{100) — CJ\f{0, 2) compound distributions using central 


difference discretisation vs the step 


r TT AT 

Size 0. Here, N — 


(I0.999/0 IS the number 


of steDS 


reouired 






5 


N 


?0.999 


time (sec) 


16 


360 


5760 


0.19 


8 


725 


5800 


0.20 


4 


1457 


5828 


0.28 


2 


2921 


5842 


0.55 


1 


5849 


5849 


1.59 


0.5 


11703 


5851.5 


5.77 


n 9^ 


23411 


5852.75 


22.47 


n 1 9^ 


46824 


5853 


89.14 


0.0625 


93649 


5853.0625 


357.03 


Table 4: 


: Example of Panjer recursion calculating the Poisson{100) — 


£A/'(0, 2; 


) compound distributions using central, forward and backward dif- 


fcreiice (liscretisatioii witli tlie step = 1. 




n 




Hn 


n 





3.72008 X 10-^^ 


2.50419 X 10-28 


1.92875 X 10-22 


1 


1.89724 X 10-^2 


5.65628 X 10-27 


2.80718 X 10-21 


5811 


0.998953196 


0.998983158 


0.998999719 


5812 


0.998953669 


0.998983612 


0.999000163 


5848 


0.9989705 


0.998999773 


0.999015958 


5849 


0.998970962 


0.999000217 


0.999016392 


5913 


0.998999942 


0.999028056 


0.999043605 


5914 


0.999000385 


0.999028482 


0.999044022 
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Table 5: Example of FFT calculating the 0.999 quantile of the 
Poisson{100) — £A/'(0, 2) compound distribution using central difference dis- 
cretisation with the step 6 = 0.5. The exact Panjer recursion for this dis- 
cretisation step gives Qo.ggg = 5851.5. 



r 


1 = 6x2'' 


0^'^ 

Vn.999 


VO.999 


Vo.999 


time (sec) 


14 


8192 


5117 


5665.5 


5851.5 


0.17 


15 


16384 


5703.5 


5834 


5851.5 


0.36 


16 


32768 


5828 


5850 


5851.5 


0.75 


17 


65536 


5848.5 


5851.5 


5851.5 


1.61 


18 


131072 


5851.5 


5851.5 


5851.5 


3.64 


19 


262144 


5851.5 


5851.5 


5851.5 


7.61 



Table 6: The weights Wi and abscissas Q of the 7-point Gaussian quadrature 

i Ci Wi 

1 -0.949107912342759 0.129484966168870 

2 -0.741531185599394 0.279705391489277 

3 -0.405845151377397 0.381830050505119 

4 0.0 0.417959183673469 

5 0.405845151377397 0.381830050505119 

6 0.741531185599394 0.279705391489277 

7 0.949107912342759 0.129484966168870 



Table 7: Convergence in DNI estimates of H[z = 5853.1) for Poisson{100)- 
£jV(0, 2) in the case of no = 1 and different truncation length K. Htaii is 
the estimate with the tail correction and H is the estimate without the tail 



correction. 


K 


H 


Htail 


time(sec) 


2 


0.9938318 


0.9999174 


0.0625 


3 


1.0093983 


0.9993260 


0.094 


4 


1.0110203 


0.9991075 


0.125 


5 


1.0080086 


0.9990135 


0.141 


10 


0.9980471 


0.9989910 


0.297 


20 


0.9990605 


0.9990002 


0.578 


40 


0.9989996 


0.9990000 


1.109 


80 


0.9990000 


0.9990000 


2.156 
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Table 8: The estimates of the 0.999 quantile, Qo.ggg, for Poisson{X)-CJ^{0, 2), 
calculated using DNI, FFT, Panjer recursion and MC methods. Standard 
errors of MC estimates are given in brackets next to the estimator. 





A 


0.1 


10 


1000 


DNI 


Qo.999 


105.36 


1,779.1 


21,149 




time 


15.6s 


6s 


25s 




K\no 


50\2 


25\1 


25\1 


MC 


\rJ r\ ooQ 
^ u.yyy 


105.45(0.26) 


1, 777(9) 


21,094(185) 




time 


3min 


3.9min 


1 1 . 7min 




Nmc 


10« 


10^ 


10^ 


Panjer 


Qo.999 


105.36 


1,779.1 


21, 149 




time 


7.6s 


8.5s 


3.6h 




h 


2-' 


2-3 


2-4 


FFT 


Qo.999 


105.36 


1,779.1 


21,149 




time 


0.17s 


0.19s 


7.9s 




h 




2-3 


2-4 




M 


214 


214 


219 



Table 9: The estimates of the 0.999 quantile, Qo.ggg, for Poisson{\)- 
GPD{1, 1), calculated using DNI, FFT, Panjer recursion and MC methods. 
Standard errors of MC estimates are given in brackets next to the estimator. 





A 


0.1 


10 


1000 


DNI 


Qo.ggg 


99.352 


10,081 


1.0128 X 10^ 




time 


21s 


29s 


52s 




K\no 


100\2 


100\2 


100\1 


MC 


Qo.ggg 


99.9(0.3) 


10, 167(89) 


1.0089(0.026) X 10^ 




time 


3.1min 


3.6min 


7.8min 




Nmc 


10« 


10^ 


10^ 


Panjer 


Qo.ggg 


99.352 


10,081 


1.0128 X 10^ 




time 


6.9s 


4.4s 


15h 




h 


2-' 


1 


1 


FFT 


Qo.ggg 


99.352 


10,081 


1.0128 X 10^ 




time 


0.13s 


0.13s 


28s 




h 


2-' 


1 


1 




M 


214 


214 


221 
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Figure 1: Panjer recursion estimate, $0,999, of the 0.999 quantile for the 
Poisson{100) — CAf {0, 2) compound distribution vs the step size 5 (top figure) 
and vs the number of steps N — $0.999 / ^ (bottom figure) . 
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Figure 2: Different approximations for the tail of the Poisson{100) — 
CM{Q, a) distribution for a) cr = 2; and b) less heavier tail a — 1. 
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